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Abstract 

We investigate two-dimensional neural fields as a model of the dy- 
namics of macroscopic activations in a cortex-like neural system. While 
the one-dimensional case has been treated comprehensively by Amari 30 
years ago, two-dimensional neural fields are much less understood. Wc de- 
rive conditions for the stability for the main classes of localized solutions 
of the neural field equation and study their behavior beyond parameter- 
controlled destabilization. We show that a slight modification of original 
model yields an equation whose stationary states are guaranteed to satisfy 
the original problem and numerically demonstrate that it admits localized 
non-circular solutions. Generically, however, only periodic spatial tessel- 
lations emerge upon destabilization of rotationally-invariant solutions. 

1 Introduction 

Neural fields (Amari 1977) describe the dynamics of distributions of activity 
on a layer of neurons. Neural fields have been suggested as models of internal 
representations in natural agents (Takeuchi and Amari 1979, Gross et al. 1998) 
as well as in robots (Steinhage 2000, lossifidis and Steinhage 2001, Erlhagen 
and Bicho 2006). Various modalities are covered such as spatial localization, 
viewing direction, attentional spotlight, the dynamics of decision making, ele- 
mentary behaviors, and positions of other agents in the environment. More ex- 
tensive studies in theoretical neuroscience (Suder et al. 1999, Mayer et al. 2002, 
Bresslofi" 2006) concentrate on primary visual cortex, cf. (Lieke et al. 1989), 

*present address: Universitat des Saarlandes, Postfach 151150, 66041 Saarbriicken 



1 



superior culliculus (Schiorwagcn and Wornor 1996), the representation of mo- 
toric primitives (Thelen et al. 2001, Erlhagen and Schoner 2002), and working 
memory in prefrontal cortex (Schutte et al. 2003, Camperi and Wang 1998). 

Recently neural fields are experiencing attention in modeling and analysis of 
brain imaging data, because they are able to represent the dynamic interaction 
of an active medium with time-varying inputs, and because the spatial and 
temporal scales in the data and neural field models are starting to become 
comparable. Especially if information about connectivity is available from the 
data (Jirsa et al. 2002) then neural fields are of high explanatory power. 

Moreover the analysis has reached a level where applications directly benefit 
from the theoretical progress, and at the same time, computational power be- 
came available that allows us to perform on-line simulations of two-dimensional 
neural fields. 

Generally speaking, neural fields serve as nonparametric representations of 
probability densities and their dynamics may perform operations on the densi- 
ties such as Bayesian computations (Herrmann et al. 1999). One-dimensional 
problems have been comprehensively studied already in the 1970s (Amari 1977, 
Kishimoto and Amari 1979, Takeuchi and Amari 1979). While for spatially 
extended, e.g. periodic patterns, the transition to the more relevant two- 
dimensional case is nontrivial but fairly well understood (Ermentrout and Cowan 
1979, Ermentrout 1998), localized activities in dimensions larger than one have 
not yet been treated with the same rigor. The situation, however, does not seem 
to be exactly complex: a large body of numerical studies together with theo- 
retical considerations (Laing et al. 2002, Laing and Troy 2003) imply a general 
instability of multi-bump solutions if the interactions are excitatory at small 
and inhibitory at large distances (see also Laing and Chow (2001) for stability 
analysis in one-dimensional models of spiking neurons). Further, there has been 
numerical evidence that single-bump solutions in two dimensions for radially 
symmetric interactions are essentially circular, which was exploited as an as- 
sumption in Taylor's early attack to the two-dimensional case (Taylor 1999). In 
ref. (Werner and Richter 2001) evidence has been provided for the existence of 
ring-shaped solutions which are possible for certain types of neural interactions. 
Along these lines one may conjecture that finite mesh-like structures of higher 
genus do exist as well. 

The situation became more spirited only recently when in Refs. (Herrmann 
et al. 2004, Bressloff 2005, Doubrovinski 2005) the stability problem of localized 
activations in two-dimensional fields was eventually tackled. Although the gen- 
erality which has been achieved in the one-dimensional neural fields is presently 
out of reach in two dimensions, a number of interesting variants of the circular 
activity configurations were analyzed so far. 

Here we present a concise and reproducible scheme for the analysis of the 
stability of localized activity distributions in neural fields. We show the ap- 
plicability of the scheme not only to the special case of circular solution but 
study also ring-shaped and bar-shaped solutions which present the complete set 
of known localized solutions for simple kernels (Werner and Richter 2001). In 
addition to the stability proofs which are based on the classical scheme (Pis- 
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men 1999, BrcsslofF 2005, Doubrovinski 2005), we are interested mainly in the 
behavior beyond the phase transitions towards the unstable regions. The desta- 
bilization of circular solution is known to lead to a transient elongation of the 
activity bubble (BrosslofF 2005, Doubrovinski 2005) which ultimately causes the 
localized solution to split or to form meandering bands. Either case is unstable 
in a strict sense: the splitting into two continues toward a plane-filling hexag- 
onal pattern while the banded patterns develop a global stripe pattern or a 
quasiperiodic arrangement. 

The destabilization is thus fundamental since for typical Mexican-hat in- 
teraction kernels there is no nearby stable state which is approached after the 
bifurcation, while the spatially extended patterns are not approached in finite 
time (unless a general criterion for convergence is drawn into consideration). Yet 
the destabilization, at least in the neurobiological applications is the most inter- 
esting part of the theory. Divergences are usually very slow and may halt com- 
pletely due to reasonable boundary conditions (cf. below), such that an activity- 
based correlational learning scheme may organize anisotropics in the connections 
which stabilizes the anisotropic activities as assumed in Ref. (Bressloff 2005) and 
exploited in (Schierwagen and Werner 1996). A theoretical account of the in- 
teraction of activity dynamics and learning was studied in (Dong and Hopfield 
1992), in relation to activity effect on feature maps cf. (Mayer et al. 2002). 

2 The neural field equation 

The neural field model describes the activations of a layer of neurons when the 
geometry of interactions rather than the specific connectivity among the neurons 
is relevant. We assume positions r G for neurons with continuous- valued 
activations u{r,t). The synaptic weights between neurons at the positions r 
and r' is expressed by isotropic interaction kernel w{r,r') = «;(|r — r'|) of 
Mexican-hat shape. Neurons are activated if their total input is greater than 
zero. We will study only equilibrium solutions without external input, and we 
neglect slow learning effects, so the synaptic weights are constant over time. 

The activation at a position results from a weighted integration over the 
inputs from all other active locations in the field and a natural decay towards 
a resting potential denoted by h. The dynamics of the neural field is thus 
determined by the equation 



where R[u] = {x \ u (r) > 0} is the excited region, i.e. a neuron receives input 
only from neurons within R. The boundary of R [u] is assumed to be smooth. 
Rescaling time, r can be set to unity without loss of generality. Equilibrium 
solutions are defined by 




R[u] 



(1) 




(2) 
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and depend on the value of the threshold parameter h and the particular form of 
w. Here we use a smooth kernel which is more general than the quasi-constant 
kernel function in Ref. (Herrmann et al. 2004). The kernel is constructed as a 
difference of Gaussian functions and is defined by four parameters K, k, M, m: 

w = Kexp(^-k\\r-r'f^ - M exp (^-m ||r - r'f ) (3) 

If non-rotationally symmetric solutions are excluded from the beginning from 
the consideration of Eq. |2] the situation simplifies dramatically and it can be 
shown (Taylor 1999) that one-bump solutions u(||r||) of certain radii are sta- 
tionary states of the dynamics ([ij. Analogously, a ring-shaped solution (Werner 
and Richter 2001) or a stripe-like solution, i.e. a degenerate ring of infinite ra- 
dius, can occur. However, when considering an arbitrarily small perturbation of 
the solution, the symmetry might be broken and new phenomena can appear, 
as will be studied in the following. 



3 Stability 

It has previously been shown that ([iJ admits rotationally invariant stationary 
solutions with disc-shaped activated region (one-bump solutions). Generically, 
these arise in the course of a "blue sky bifurcation" (Strogatz, 1994): no solution 
is present in the subcritical parameter region whereas two solution branches 
bifurcate as control parameter exceeds the critical value. The two solutions 
are rotationally invariant one-bumps. Stability analysis of these states with 
respect to rotationally invariant perturbation is essentially equivalent to stability 
analysis of one-bump solution of the one-dimensional model. It reveals that the 
unstable branch generically corresponds to the bump of smaller radius. Upon 
destabilization the region of activation expands as the solution approaches the 
stable branch, corresponding to the circular bump of larger radius (See Fig. 111. 
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Figure 1: (a) Stable one-bump solution of the neural field equation. Parameters 
are K = 2.5, k = 5, M = 0.5, m = 0.5, h = -0.281 (b) Unstable one-bump 
solution for the parameter values K = 2.5, k = 5, M — 0.5, m — 0.5, h = —0.294 
(c) Annular solution with the parameters K = 2.5, k = 5, M = 0.5, m = 0.5, 
h — —0.115 (d) A part of a stable stripe-like solution. Parameters are the same 
as in (c). At the same parameters also an unstable solution exists (not shown). 



Seeking a stationary solution u of the two-dimensional model ([ij assuming 
rotational invariance of the field (i.e. il(r) = u(r) with r — ||r||) leads to a 
problem, essentially equivalent to that of finding stationary states in the one- 
dimensional model (see Appendix). Two types of solutions besides the circular 
one-bumps are readily constructed: solutions with annular activated regions and 
solutions with stripe-shaped region of activation (see Fig. [iJ. The possibility of 
existence of the former has been pointed out previously, whilst the latter can 
be seen as degenerate annuli of infinite inner radius. 

We now turn to the analysis of the stability properties of the above- 
mentioned stationary states. Consider the one-bump solution with disc-shaped 
activated region (i.e. u{r) > iff r < R ). Consider dynamics of a small 
perturbation trj, i.e. 

u{v,t)=u{\\v\\) + tr^{r,t) (4) 

Inserting into ([T]| and keeping terms of order at most 1 in e one finds that to 
linear order the dynamics of the perturbation obeys 



dt 



W I 



\)5{u{y'))r]{v',t)dv' 



(5) 
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where 5 (u(r')) is Dirac delta function (see Appendix). Substituting an Ansatz 
of the form ij (r, t) = exp (Xt) ^ (r) one arrives at an eigenvalue problem which 
in polar coordinates becomes (see Appendix) 

Ae (r, 6) = -e (r, 0) + r / g {9 - 9') ^ (r, 9') d9' (6) 



Here g is a 27r periodic function depending on kernel t(;(||r — r'||) and F is a 
constant given by 



T = R 



du (r) 




dr 


r=R 



(7) 



i.e. the ratio of the radius of activated region R to the absolute value of the 
slope of the radial profile of the stationary solution at r = i?. Clearly, explicit 
evaluation of F requires calculating the stationary solution which is implicitly 
given by ^ in terms of a double integral. Equation (|6| is known as Fredholm's 
integral equation of the second kind. The integral operator in the right-hand side 
of ([6} is compact, bounded and self-adjoint implying that every spectral value 
is an eigenvalue, all eigenvalues are real, each eigenspace is finite-dimensional 
and zero is the only possible accumulation point of eigenvalues (Kreyszig, 1978). 
Solving ([6} we obtain eigenvalues and eigenfunctions which in polar coordinates 
read (see Appendix) 

A„ -1 + F /q g (9) cos in9) d9 

e„ = Ja'' w (r, 9, R, 9') cos {n9') d9' (8) 

where R is the radius of the circular activated region. The nth eigenfunction 
is 27rn-periodic in 6', implying that it is D„-symmetric (symmetry with respect 
to rotation by 27r/n around the origin and with respect to reflections on the n 
respective symmetry planes; the shape in Fig. [5^, e.g., is D4-symmetric) and 
corresponds to a multi-periodic deformation of a circle. Eigenvalue spectra for 
one-bump solutions are given in Fig. [2] Certain features of these are readily 
interpretable. E.g. we see that for the bump of smaller radius the eigenvalue Aq, 
corresponding to rotationally invariant deformation, is positive implying insta- 
bility with respect to perturbations in radius of the bump. The spectrum, cor- 
responding to the larger bump, however, is non-positive, implying stability with 
respect to arbitrary perturbation in agreement with earlier results. Also, the 
eigenvalue Ai, that corresponds to a 27r-periodic deformation (or, equivalently, 
to a translation of the bump) vanishes, reflecting the translational invariance of 
([ij. Exploiting this observation, it appears possible to re-express F in Q more 
explicitly as F = 1/ J^^ g (9) cos (9) d9 , whereby the expressions for the other 
eigenvalues simplifies to 

,^_,j:i9i0)cosin9)d9 
J^" g (9) cos (9) d9 
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Figure 2: Spectra of corresponding to solutions in Fig. [TJa, b, and c, respec- 
tively. The positivity of the zeroth mode in (b) indicates the instability of this 
solution with respect to circular perturbations. The solution (a) is stable, with 
a marginal instability with respect lateral shift as visible from the vanishing first 
eigenvalue. The instability of the annulus solution in (c) shows a characteristic 
scale which causes the solution to split into three single bumps. 



Note that contrary to ([8|, Eq. |9]does not contain explicitly the stationary 
solution M, which allows us to calculate the nth eigenvalue as a function of the 
radius of activated region without evaluating double integrals in the implicit 
expression for u. Only integrals over one-dimensional manifolds appear in ([o}, 
greatly simplifying the calculation of the spectrum. Apart from theoretical con- 
siderations, this is of importance for technical applications since the knowledge 
of stability properties of solutions of ([ij could affect their use for representing 
probability distributions e.g. in implementations of autonomous robot memory. 




a R b 



Figure 3: (a) Four eigenvalues versus bump radius R. The curve which repre- 
sents Aq described the stability with respect to perturbations in bump radius. 
The others are (from left to right) A2 (reflection-invariant deformation), A3 (-D3- 
invariant eigenmode) and A4 (D4-invariant eigenmode). Ai is identically zero 
reflecting the metastability of the solutions w.r.t. lateral shifts. Parameters are: 
= 1.5, fc = 5, M = 0.5, m = 1.5 (b) Spectrum determining stability of 
stripe-shaped solutions shown in Fig. [l|d). Only the information for the stable 
solution if given here. This solution looses stability for a band of values of Q.. 

Previous work (Werner-Richter 2001) conjectured the bifurcation branch cor- 



7 



responding to the stable one-bump to remain stable for all values of control 
parameter. Using (|9} this is readily proved false: higher and higher frequency 
eigenmodes progressively turn unstable as the radius of the stationary solution 
is increased (see Fig. [3]). 

Strictly speaking, the assertion that linear stability analysis as outlined 
above correctly determines stability properties of stationary solutions relies on 
additional assumptions on operators appearing in the right-hand side of ^ 
(Schaefer and Golubitsky 1988). For infinitely-dimensional non-smooth fields 
these generically need not hold. In order to check whether stability analysis 
is indeed adequate, correctly determining stability properties of the stationary 
solutions, we performed a number of numerical simulations. Fig. Q depicts a 
simulation of unstable one-bump whose corresponding stability spectrum reveals 
that the maximal (positive) eigenvalue is that of the Z32-symmetric eigenmode. 
As the time elapses, initially circular activated region keeps deforming, forming 
blob-like protrusions. These subsequently bud off from the middle-bump. The 
newly-formed activated domains keep splitting, progressively tiling the plane in 
a hexagonal pattern. We would like to stress that the pattern which is formed 
immediately after destabilization of the stationary state is I?2-symmetric as 
would be expected from the properties of the eigenvalue spectrum. That is, sta- 
bility properties of the stationary solutions as well as qualitative aspects of the 
pattern, forming upon destabilization of the steady state, are correctly predicted 
from the above-mentioned linear stability analysis. 






"9? BB 




Figure 4: (Upper panel) Time evolution of an unstable bump, undergoing 
stability loss at a _D2-eigenmode. Parameters are: K = 1.5, fc = 5, M = 0.5, 
m = 1.5, h = —1.46 • 10~^. (Lower panel) Time evolution of an unstable bump, 
undergoing stability loss at a invariant eigenmode. Parameters: K = 1.5, 
k = 5, M = 0.5, m = 1.5, h = -4.43 • lO^^. 



Symmetry breaking accompanying destabilization of a stationary bump was 
examined for a broad range of parameters (e.g. parameters which yielded solu- 
tions with maximal eigenvalue of the respective spectrum corresponding to D3-, 
D4- and Dg-symmetric perturbations, see Fig. |4]|. In all of the cases, the course 
of the symmetry breaking was appropriately determined by linear stability. 

Stability analysis of annular solutions proceeds along the same lines as that 
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of one-bumps (see Appendix). The essential difference is that instead of a singe 
equation ([g}, a system of two equations results, meaning that to every non- 
negative whole number corresponds a pair of real eigenvalues. (In general, to 
every boundary of an activated domain there corresponds an equation in the 
corresponding eigenvalue problem). We find that in the case of annular solu- 
tion shown in Fig. [ijl the largest eigenvalue corresponds to a Z?3-symmetric 
perturbation. Simulations demonstrate that initially rotationally invariant an- 
nulus splits into three adjacent blobs which gradually drift apart (not shown). 
The symmetry of the resulting state is the same as that of the largest eigen- 
value. Again, stability of the stationary solution as well as qualitative aspects 
of the emerging pattern are readily predicted from the analytically computed 
spectrum. 

Finally, let us consider the stripe-shaped solutions. Their stability is gov- 
erned by an eigenvalue problem, similar to that governing stability of the annuli. 
However, the corresponding linear operator is no longer compact (this is a conse- 
quence of activated region being unbounded) and the spectrum needs no longer 
remain discrete. Actually, in this case the spectrum is continuous: to every 
real corresponds a pair of (real) eigenvalues. Fig. [sja shows results of linear 
stability analysis of the stripe solution, depicted in Figure p. In the corre- 
sponding simulation the stripe is seen to split into a row of separate bumps - 
a "chain-of-pearls" configuration. The inter-bump separation is the same as the 
wave-length of the eigenmode, corresponding to the largest eigenvalue. Again, 
stability and semi-quantitative properties of patterns resulting from the station- 
ary state destabilization are readily predicted from the respective eigenvalue 
spectrum. 

In summary, preceding section describes all of the stationary non- 
homogeneous solutions of the two-dimensional Amari equation known up to date 
and exhaustively examines their stability properties. Quite strikingly, stability 
analysis of the non-homogeneous steady state is possible since the eigenvalue 
problem ([sf is effectively one-dimensional although a two-dimensional system is 
being considered. 
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Figure 5: Dynamics of the modified field with and without noise, (a) In the 
absence of noise an unstable bump develops four identical protrusions subse- 
quently splitting into four separate blobs, (b) In presence of noise a ring-shaped 
region of activation is initially annular, but (c and d) develops gradually into 
an irregularly meandering profile. Parameters are k — 5, m = 1.6, K = 1.5, 
M = 0.45, h = -0.0505, p = -6.25 • 10"^ 



4 Modified equation 

A long-standing question regarding Amari-model is existence of non-rotationally 
invariant stationary solutions with bounded and connected region of activation. 
Such states are likely to bifurcate from circular solutions upon destabilization 
at a non-rotationally invariant eigenvalue. In fact, stability loss of one branch 
is always accompanied by emergence of another in its vicinity provided that the 
mappings defining the dynamical system under consideration are sufficiently 
smooth (Crandall and Rabinowitz 1971). However, in the above simulations 
exclusively periodic patterns resulted upon destabilization of circular solutions. 

In order to address existence of non-rotationally invariant solutions of ^ 
with bounded and connected activated region we shall modify the original 
Amari-model ([T]| so as to obtain a related ("modified") equation fulfilling the 
following three conditions. 

1. Stationary solutions of the modified equations should be solutions of the 
unmodified equation 

2. The modified equation should not admit solutions with unbounded acti- 
vated region. 
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3. Stability of a solution of Eq. ([2} should remain unaltered by the modifi- 
cation. 



According to condition 3 rotationally invariant stationary solutions behave like 
those of the unmodified Amari-model, admitting symmetry breaking at a non- 
rotationally invariant eigenvalue. However, the symmetry breaking cannot re- 
sult in a spatially extended periodic pattern according to condition 2. Conse- 
quently, a non-rotationally invariant localized state is likely to emerge. It will 
be a solution of the original (unmodified) Amari model with desired properties, 
provided that its region of activation remains connected in the course of desta- 
bilization. Note, however, that conditions 1-3 do not suffice to ensure that the 
activated domain will not start splitting into separate disconnected region upon 
symmetry breaking. 

We now turn to the construction of a modification of ([ij satisfying the above- 
mentioned conditions 1-3. Consider some circular one-bump solution Uh of (jlj 
with resting potential h and the area of activated region ^[w^,]. Let us modify 
([iJ according to 

dtu=-u{r,t)+ [ [w{\r-r'\) + q]dr' -qA[uh] + h' (10) 

Jr[u] 

where q is any real number. Suppose that u is a stationary solution of ([To| for 



a certain q. Substituting u into (10 1 one obtains with the area of R[u] being 
denoted by A[u]: 



= -u{r,t) + [w{\r-r'\) + q] dv' - qA[uh\ + h' ^ 
= -u (r, t) + w (|r - r'l) dv' + qA[u\ - qA[uh] + h' 



(11) 



implying that m is a solution of Q with the resting potential h replaced by 
h' — qA\uh] + (7A[m]. Consequently, condition 1 is satisfied. Note that the 
circular one-bump solution Uh of the original problem ([2} which was used when 
constructing ([To| solves the modified problem ( 10 1 with h' — h and any q. 



Eq. ( 10 1 does not admit stationary solutions with unbounded region of ac- 
tivation. Indeed, assuming that such a solution u exists, substituting ii into 



(101 we were to conclude that the integral term of the right-hand-side of ( 10 1 is 
infinite which forms a contradiction. Therefore condition 2 above is satisfied. 



Finally, we show that the modification (10 1 preserves stability properties 



of the stationary solution (which is a stationary solution of both the mod- 
ified and the unmodified problems ([Tl and (101 respectively by construction). 
Recall that in deriving equation (|8| we did not make use of any particular 
assumptions on the form of the integral kernel i(;(|r — r'|). Consequently, 
this expression for the eigenvalue spectrum is equally valid for the modified 
problem as well as for the unmodified one. Note, however, that when de- 
riving stability spectrum in the case of the modified problem ( [To| we shall 
exchange g {6 — 0') by g{9 — 0') + qA\uh\ (see derivations in the Appendix). 
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Using J^^ qA[uh] cos {n9) dO — qA\uh\ J^^ cos {n9) dO — gA[?I;j]27r(5„o it now fol- 
lows from ([sj that all eigenvalues of the stability spectrum of Uh (except for 
Ao corresponding to perturbations of the radius of the bump) remain unaltered 
by the modification. Consequently, if Uh is unstable with respect to some non 
rotationally-invariant perturbation in the original problem dTl, it is unstable 
with respect to such a perturbation in the modified problem (10 1, whereby con- 
dition 3 holds. 

The above arguments imply that a rotationally invariant solution Uh of 
Eq. (jlj that is unstable at a non-circular eigenmode solves also the modified 
equation ( 10 1 and is unstable with respect to the same eigenmode of the dy- 
namics (|10|. Furthermore, contrary to the case of ([ij, the dynamics of (10 1 can 
never result in a periodic pattern with unbounded activated region if g < 0. 

As stated above, conditions 1 - 3 do not suffice to guarantee that the acti- 
vated region will remain connected as ut follows the dynamics (10 1. Neverthe- 
less, by tuning the parameter g in (10 1 one is able to trap the dynamics in the 
vicinity of instability in a state with connected activated region. 





Figure 6: (a) Contour plot of the I?2-symmetric solution superimposed with 
its level-curves. The zeroes level-curve (corresponding to u = 0) delimits the 
activated region. The symmetries of the equation and initial condition allow 
to restrict the simulation to one quarter of the domain. The resulting blob is 
roughly ellipse-shaped with the vertical semi-axis being longer than the horizon- 
tal one. (b) Cross-section through the profile of the stationary solution along 
the coordinate axes, (c) Rate of growth of the solution versus iteration step cal- 
culated as maximal deviation between two subsequent time steps. Parameters 
are K = 1.5, k = 5, M = 0.5, m = 1.5, h = 1.50, q = 0.02850. 



Assume that for q = the dynamics in vicinity of the bifurcation tends to 
increase the area of the activated region. Note that q could be understood as 
a Lagrange multiplier, which ensures the constancy of the area for one special 
value, which we certainly exclude. Suppose now that q is chosen to be negative. 
The additional term in the integrand of ( 10 1 will tend to counterbalance the 
area increase. We can now choose q such that two effects counterbalance and 
a non-rotationally invariant steady state with bounded and connected region of 
activation will result. The calculation of effective values of q requires the consid- 
eration of higher orders of the dynamics which can be expressed by a Ginzburg- 
Landau equation, cf. Doubrovinski (2005). Here, we will instead take resort to 
numerical simulations. These confirm the emergence of non-rotationally invari- 
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ant steady states with bounded and connected region of activation. For example, 
Fig. |6] shows the destabiHzation of a circular one-bump which is unstable at a 
Z?2-symmetric eigenmode, developing into a non-rotationally invariant steady 
state with ellipse-shaped activated region. Only one quarter of the domain was 
simulated (the field in the other three quadrants is determined by that on the 
simulated quadrant due to Euclidean symmetry of dynamic equations and D2- 
symmetry of initial conditions) on a grid of 300 x 300 pixels in order to increase 
the accuracy. Symmetries corresponding to higher eigenvalues are irrelevant 
because the eigenvalues A„ are stable for n > 3 for the given parameters. The 
length difference of the half-axes of the activated domain of the resulting D2- 
symmetric stationary state was much larger than the spatial discretization step 
(some 30 pixels) allowing to conclude that that the stationary solution obtained 
is not a discretization artifact. 

Another example is given in Fig. [7] showing the dynamics in the vicinity of 
stability loss at a I^a-symmetric eigenmode. Initially the field is Da-symmetric. 
We see that initially a Da-symmetric state with connected activated region does 
indeed result. The system dwells in this state during considerable time, breaking 
Da-symmetry due to small numerical perturbations (in absence of perturbations 
Da-symmetry must be preserved by the dynamics (10 1 due to invariance under 
Euclidean transformations). Supposedly, the emerging Da-symmetric state is 
stable when the dynamical system is confined to the subspace of Da-symmetric 
functions. The final peak in Fig. [7^, however, indicates an instability with 
respect to a D2-symmetric perturbation. 



□ □ 






Figure 7: (a)-(d): Time evolution of activated region, (e) Rate of growth of 
the solution versus iteration step calculated as maximal deviation between two 
subsequent time steps. The initial deep "valley" corresponds to the shape of 
activated region shown in panel (b). Parameters are: K = 1.5, k = 5, M = 0.5, 
m = 1.5, h = -0.0260, q = 0.0125. 



In conclusion we would like to stress that the particular choice of modifica- 



tion according to (10 1 is very restrictive. In fact, many other equations whose 
stationary states satisfy ([ij and do not admit solutions with unbounded re- 
gions of activation are readily constructed along the same lines. For instance 



dtu = -u + /j 



r — r'l) dr' — /^r^i |r — r'| dr'dr + h can be shown to have 
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these properties, arguing essentially as when proving that conditions 1 and 2 
above are satisfied by stationary states of (10 1. We believe that further investi- 
gation of such modifications will provide insights into the properties of unstable 
solutions of Amari-equation. 



5 The time course of symmetry breaking 

Bifurcation theory describes the time course to critical behavior in low- 
dimensional systems. Under certain conditions the parametric destabilization 
of an activity distribution does not lead to a nearby stable state, but initiates a 
cascade of symmetry breaking events which eventually approaches a new distant 
stable configuration. It has been a motivation for the present study to demon- 
strate the complex evolution of the state of the field after the loss of stability. 
It is these configurations that bear the greatest computational potentials. 

For the Amari equation with a simple kernel only the existence of rotation- 
ally invariant bumps has been proven. This work demonstrates existence of two 
other non-periodic stationary states: stripe-shaped and annular solutions. Sta- 
bility analysis of these would be expected to be very involved. Surprisingly, the 
special form of the Amari equation makes this stability problem amenable to an- 
alytical treatment. The total synaptic input to a given neuron is only dependent 
on the shape of the boundary of the activated domain, making the eigenvalue 
problem effectively one-dimensional, thereby allowing for fairly straight-forward 
calculation of the spectrum for each of these cases. 

Strictly speaking, spectral properties of the linearized operator do not guar- 
antee the stability of the stationary state unless additional assumptions are 
satisfied. However, our extensive numerical investigation shows that stability is 
indeed correctly predicted by eigenvalue analysis. Furthermore, the spectrum 
allows to predict certain semiquantitative features of solutions approached after 
the onset of instability. 

Our numerical experiments showed that exclusively spatially extended so- 
lutions (i.e. those with unbounded activated domain) appeared in the unstable 
parameter region. Yet, a slight modification of the interaction kernel, introduc- 
ing long-range interactions circumvents this, yielding non-rotationally invariant 
stationary states with bounded and connected region of activation which at the 
same time are stationary states of the original unmodified equation. This set- 
tles a longstanding question regarding existence of solutions of this type. This 
bears also relevancy for biological systems. Assuming a certain degree of shift- 
twist symmetry the existence of elongated blobs suggests mechanisms for the 
emergence of orientation selectivity in neurons of the primary visual cortex. 
Although the degree of asymmetry of non-rotationally invariant solutions was 
moderate, these effects could in principle be enhanced by (Hebbian) learning 
which we disregarded in our treatment. 
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A Appendix 

A,l Circular one-bump solution 

Consider the development of a small perturbation erj of a stationary rotationally- 
invariant one-bump u. Substituting into the Amari equation and linearizing in e 
we have 



dy (x, t) 
dt 



= -?7(x,i)+ / u;(|x-x'|)(5(zl(x'))77(x',i)dx'. (12) 

JR2 



In polar coordinates we write (somewhat informally) ti;(|x — x'|) 
w (r, 9, r' , 9'), and use the rotational invariance of u, such that 



dv{r,9,t) '■'^ 
-^^ = -v{r,9) + 

Recall that 



"'0 



r'w (r, 9, r', 9') S {u (r')) V (r', 9') dr'd9' (13) 



5 


{x - Xi) 










dx 





(14) 



where the sum is over the roots of /, provided that / is differentiable at the 
corresponding points. Using ([T4|, (|2l simplifies to 



dij (r, 9) 
dt 



27T 



^~ri{r,t)+ / Rw{r,9,R,9') 



1 



du{R) 
dr 



r]{R,9')d9', (15) 



where R is the radius of the (circular) region of activation of u. Substituting 
r] = e^*^ (r, 9) we arrive at the following eigenvalue problem. 



27r 



(r, 9) = (r, 0) + r / w (r, 9, R, 9') ^ [R, 9') d9' , 



(16) 



where T = R/ \{du{R) /dr)\. By restriction of both sides to r = i?, the eigen- 
value problem (IT6| is solved by 



A. 



^„ (i?, 9) = cos n9 
4 + r Q"' w {R, 9, R, 0) cos {n9) d9, 



(17) 



where n are nonnegative integers. The last equation can be verified by noting 
that w {R,9, R,9') is a function of {9 — 9') alone and Fourier-expanding the 
integrand of (161. 

The r-dependence of the eigenfunctions can be derived from (17 1 by exploit- 
1 of the 

^n{r,9) 



ing the special form of the eigenvalue problem (16 1. 



w{r,9,R, 9') COS {n9')d9' 



(18) 
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Shift invariance allows us to conclude that Ai = 0. Thus, T is obtained more 
explicitly from Eq. W7\ 



r = 



1 



/o w (i?, R, 0) cos {9) dO 



(19) 



Now the eigenvalues can be calculated without evaluating the slope of the sta- 
tionary solution u: 



/p w {R, 9, R, 0) cos {n9) d9 
J^"" w{R,0,R,O)cos{0) d9 



(20) 



Eq. (20 1 is particularly convenient for examining stability properties of circular 



one-bump solutions. 



A,2 Annular solutions 

The existence of solutions with annular region of activation was suggested al- 
ready in (Amari 1977). Denoting the inner radius by i?i and the outer radius 
by i?2, ^15]) is immediately rewritten 



-77 (r, 9) + Ti w (r, 0, Ri,9') 0') d9' 

-T2 / w{r,e,R2,9')7j{R2,e')d9', (21) 
Jo 



drj (r, 9) 
dt 



where Ti ^ RJ \du (Ri) /dr\, T2 = R2/ \du {R2) /dr\. Analogously to the 
derivation of (161, we set 77 = e'*'*^ ('^i^') a-nd restrict both sides to Ri and i?2- 



Aa = -Ci + Ti Jo " w {Ru9, Ri,9') ^^9' + T2 " w {Ri,9, R2, 9') C^d^' 
= -^2 + Ti /o " w {R2, 9, R,,9') Cid9' + T2 " w {R2,9, R2, 9') ^^d0'(22) 

where = £^{Ri,9) and — £,i{9'), i E {1,2}. It is natural to seek the 
eigenfunctions of the form [^i,'C2] = vcos {n9), where v is some two-dimensional 



vector. Using this substitution it turns out that the eigenvalues of (22 I are the 
same as those of the matrices 



-1 + Fi J^^ wii cos n9d9 T2 J^^ W12 cos n9d6 



Ti J^^ W12 cos n9d9 — 1 + r2 J^^ W22 cos n9 



2-iT 



such that Eq. 23 allows us to construct the spectrum of (22 I. 



(23) 



A. 3 Stripe-shaped solutions 

Stripe-shaped solutions can be constructed by assuming the region of activation 
of the form R[u] — {{x,y) \ Q <x < L} and can be seen as degenerate annuli 
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with infinite inner radius, considered in the previous section. Interestingly, in 
this case it appears possible to give an explicit expression for the stationary 
one-bump solution in terms of the model parameters: 



(24) 



u(x,y) = 2^{Km eviy^kxj — kM erf (-y/rria;) — 
-Km erf (^-^/kL+Vk xj + Mk eri{—y/mL + y/mx)) 



The counterpart of (21 1 now reads 



dr] 

di 



1]+ / w{\x~x'\)5{u{x'))f]{dx') 



w {x, y, 0, y') ^^^^^ rj (0, y') dy' 

dx 



w{x,y,L,y') 



1 



du{L) 
dx 



v{L,y')dy', 



(25) 
(26) 



where (a;,?;) and {x',y') are Cartesian coordinates of x and x' respectively. As 
in the former cases, substituting ri{x,y,t) = e^*^ (x, y) and restricting to or 
to L we find 



ACi = -6 + Ti JZ, w (0, 0, 0, y') Cidy' + JZ, ^ (0, 0, L, y') ^',dy' 

A6 = -6 + Ti JZ, w (L, 0, 0, y') i[dy' + fZ, w (0, 0, L, y') ^',dy' (27) 



where now denotes (y'). Subscripts 1 and 2 designate restrictions to a:; = 
and X = L, respectively, according to Fi = 1/ \du{0) /dx\, T2 = 1/ \du{L) /dx\ 
and ^1 (y) = ^(0,y), ^2(2/) = ^{L,y)- Arguing in exactly the same way as 
when considering annular solutions we conclude that ([26} admits an uncountable 
infinity of eigenvalues which are the same as those of matrices 



T2WL m 



(28) 



where woi^) = /^^ (0, 0, 0, y') cos (f^y') rfy' 

JZ^w{L,0,L,y')cos{ny')dy', = JZ^w{L,0,0,y')cos{ny') dy' = 

JZ w (O7 L, y') cos (riy') dy'. To every nonnegative real number corresponds 
a pair of eigenvalues. The corresponding eigenfunctions can be evaluated 
from (|27l). 
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